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Abstract 

^ It is a commonly observed phenomenon that spherical particles with inertia in an 

incompressible fluid do not behave as ideal tracers. Due to the inertia of the parti- 
cle, the dynamics are described in a four dimensional phase space and thus can differ 
considerably from the ideal tracer dynamics. Using finite time Lyapunov exponents 
Q we compute the sensitivity of the final position of a particle with respect to its initial 

^) velocity, relative to the fluid and thus partition the relative velocity subspace at each 

point in configuration space. The computations are done at every point in the relative 
velocity subspace, thus giving a sensitivity field. The Stokes number being a measure 
of the independence of the particle from the underlying fluid flow, acts as a parameter 
in determining the variation in these partitions. We demonstrate how this partition 
framework can be used to segregate particles by Stokes number in a fluid. The fluid 
model used for demonstration is a two dimensional cellular flow. 
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oo 1 Introduction 

o 

It has long been observed that particles with a finite size and mass have different dynamics 
from the ambient fluid. Because of their inertia the particles do not evolve as point-like 
tracers in a fluid. This leads to preferential concentration, clustering and separation of par- 
ticles as observed in numerous studies [H El E]. The inertial dynamics of solid particles 
can have important implications in natural phenomena, e.g., the transport of pollutants 
and pathogenic spores in the atmosphere [U |5], formation of rain clouds [6] by coalescence 
around dust particles and formation of plankton colonies in oceans [7j. Similarly, the inertial 
dynamics of reactant particles is important in the reaction kinetics and distribution of reac- 
tants in solution for coalescence type reactions [8] . Mixing sensitive reactions in the wake of 
bubbles has been shown to be driven by buoyancy effects of reactants [9]. Recently, a prin- 
ciple of asymmetric bifurcation of laminar flows was applied to the separation of particles 
by size and demonstrated the separation of flexible biological particles and the fractional 
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distillation of blood [TQl [TT] . Innovative channel geometries have been empirically designed 
to focus randomly ordered inertial particles in microchannels [12]. These phenomena and 
related applications rely on the non-trivial dynamics of inertial particles in a fluid. In this 
paper, we demonstrate a theoretical tool to achieve particle segregation, by studying the 
sensitivity of the dynamics of inertial particles in a fluid. 

We employ a simplified form of the Maxey-Riley equation [15] as the governing equation 
for the motion of inertial particles in a fluid. The dynamics of a single particle occur in 
a four dimensional phase space. The sensitive dependence of the particle motion on initial 
conditions is quantified using the finite time Lyapunov exponents (FTLE). It has been shown 
previously [13l E] , that the ridges in the FTLE field act as separatrices. These are in general 
time dependent and go by the name of Lagrangian coherent structures (LCS). We chose to do 
a simplified sensitivity analysis by perturbing the initial conditions in only two dimensions, 
in the initial relative velocity subspace. We obtain a sensitivity field akin to a FTLE field but 
restricted to the relative velocity subspace and demonstrate numerically that the ridges in 
this field act as separatrices. The partitions in the relative velocity subspace created by these 
separatrices determine the eventual spatial distribution of particles in the fluid. Using this 
partitioning scheme we show how the Stokes number acts as a parameter in the separation 
of particles of different inertia or size. 

The paper is organized as follows. In §2 we review the the equation governing the 
inertial particle dynamics in a fluid and it's simplified form. In §3 we briefly review the 
background theory of phase space distributions of finite time Lyapunov exponents, which we 
use to quantify the sensitivity of the physical location of inertial particles with respect to 
perturbations in the initial relative velocity. We also describe our computational scheme to 
obtain the sensitivity field in the relative velocity subspace. In §4 we present results for the 
sensitivity field of the inertial particles in a cellular flow. In §5 we demonstrate our procedure 
for the segregation of particles by their Stokes number using the results from §4. In §6 we 
give numerical justification for the robustness of the sensitivity field to perturbations in the 
velocity field of the fluid. In §7 we discuss the results and give conclusions. 



2 Governing Equations 

Our starting point is Maxey and Riley's equation of motion of a rigid spherical particle in a 
fluid ([15]). 

PpTt = Pt ^ + (*p-P/)g-^?(v- u--V 2 u) (1) 

- P/ (---(u--V 2 u)) 
Hn dt Dt K 6 

9p/ fv f* 1 d , a 2 ~ . , 
—\ - / — (v-u V 2 uWr 

where v is the velocity of the solid spherical particle, u the velocity field of the fluid, p p 
the density of the particle, pj, the density of the fluid, v the kinematic of the viscosity of 
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the fluid, a, the radius of the particle and g the acceleration due to gravity. The term on 
the right hand side are the force exerted by the undisturbed flow on the particle, the force 
of buoyancy, the Stokes drag, the added mass correction and the Basset-Boussinesq history 
force respectively. Eq ([!]) is valid under the following restrictions. 

a(v-u)/L « 1 (2) 
a/L « 1 

<7><I> « 1 



where L and U/L are the length scale and velocity gradient scale for the undisturbed fluid 
flow. The derivative 

Du du , „ s 

^ = ¥ + KV)u (3) 
is the acceleration of a fluid particle along the fluid trajectory whereas the derivative 

<^ v <9v . _. . t . 

is the acceleration of a solid particle along the solid particle trajectory. 

Eq. can be simplified by neglecting the Faxen correction and the Basset-Boussinesq 
terms [T7]. We restrict our study to the case of neutrally buoyant particles, i.e p p — pf . 
Writing W = (v - u), the relative velocity of the particle and the surrounding fluid, the 
evolution of W becomes 

^ = -(J + /i/)-W (5) 
and the change in the particle position is given by 

^ = W + u (6) 

where J is the gradient of the undisturbed velocity field of the fluid, u, and p, = ^St^ 1 is a 
constant for a particle with a given Stokes number St. Eqs (|5| and ^ can be rewritten as 
the vector field 

with x = (r, W) = (x,y,W x ,W y ) G M 4 . Eq (J7h defines a dissipative system with constant 
divergence — It has been shown by Haller [18] that an exponentially attracting slow 
manifold exists for general unsteady inertial particle motion as long as the particle Stokes 
number is small enough. For neutrally buoyant particles this attractor is W = (the xy 
plane). Despite the global attractivity of the slow manifold, domains of instability exist in 
which particle trajectories diverge [HI [161 HH]- 
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3 Sensitivity Analysis 



The Lyapunov characteristic exponent is widely used to quantify the sensitivity to initial 
conditions. A positive Lyapunov exponent is a good indicator of chaotic behavior. We have 
used the finite time version of the Lyapunov exponents, the FTLE, a measure of the the 
maximum stretching for a pair of phase points. 

We review some important background regarding the FTLE below, following [T31 Q3J . 
The solution to eq ^ can be given by a flow map, 0^ o , which maps an initial point x at 
time to t° x t at time t. 

x* = <(x ) (8) 

The evolution over a time T of the displacement between two initially close phase points, 
x(t ) and x(t ) + 8x.(t ), is given by 

5x(t + T) = # ^ + J (x) (5x(t ) + 0(||<$xf) (9) 
Neglecting the higher order terms, the magnitude of the perturbation is 

\Mto + T)\\ 



\ 



^(x)*^(xj 



M* ), \ rfx * x(to) ' (10) 



The matrix 



c #r T (x) ^ < +r (x) 

dx dx 

is the right Cauchy Green deformation tensor. Maximum stretching occurs when the per- 
turbation <5x is along the eigenvector n max corresponding to the maximum eigenvalue A max 
of C. The growth ratio is given by 

||5x(to + T)||/||5x(t )|| =e ffl ^»l T l (12) 

where 

<Tl(x(t )) = p In y/XmaxiC) (13) 

is the maximal finite time Lyapunov exponent. One can associate an entire spectrum of 
finite time Lyapunov exponents with x(£o), ordering them as 

<7i(x(t )) > <r 2 (x(t )) > CT 3(x(t )) > ^(x(t )) (14) 

The entire spectrum of the Lyapunov exponents can be computed from the state transition 
matrix using singular value decomposition. 

$(t, t ) = B{t, t )A(t, t y/ 2 R(t, t ) (15) 
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The diagonal matrix A gives all the Lyapunov exponents. 



S(t / ,t ) = ln[A(t / ,to) 1/im ] (16) 

where T = tf — t and £(t/,to) = diag(ai, ...,(74). An arbitrary perturbation in the fixed 
basis can be transformed using a time dependent transformation 



5x'{t) = A(t,t ,tf)6x(t) (17) 
such that in the new basis (the primed frame), the variational equations become 

5±'(t) = E(t f ,t )5x'(t) (18) 
Since H(tf,to) is a constant diagonal matrix, we have 

5x'(t) = e (t -' o)s( '" io W(t ) (19) 

The first coordinate in the new frame grows as 5x[(t) = e^ t ~ to > (T1 5x' x {to) . The time 



dependent transformation A(t) is given by [2Tj . 

A(t,t ,t f ) = e^^f^Ritf^toTRit^oMt^or^Bit,^) (20) 



3.1 Sensitivity to initial relative velocity 

Since the dynamics of the inertial particle is in a four-dimensional phase space, the separa- 
trices, that is LCS defined by ridges in the field of the maximal FTLE, are three dimensional 
surfaces (see [H]). However, because the system is dissipative and the global attractor is 
the xy subspace, we can obtain meaningful information by restricting the computations to 
a lower dimensional subdomain of the phase space. This we do by considering an initial 
perturbation only in the relative velocity subspace and study how this perturbation grows 
in the xy plane, the configuration space, i.e, 

5x(t ) = [0,0, AW x ,AW y }* (21) 

where AW X , AW y are the perturbations in the relative velocity subspace. Using the time 
dependent transformation A(t,t ,tf) the evolution of the perturbation is given by 

5x(t) = A~ 1 (t)e (t -* o)s ^'* o) A(to)5x(to). (22) 

The growth of perturbation in the xy plane is given by the first two components of 
the above vector. The last two components of the above vector are the evolution of the 
perturbations in the relative velocity subspace. Since the xy plane is a global attractor 
these tend to zero. One can choose a finite time, T, such that the evolution of the initial 
perturbation comes arbitrarily close to the xy plane. In this way the sensitivity of the final 
spatial location of the particles to initial relative velocity can be computed. 
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3.2 Numerical computation of the Sensitvity Field 

As equation(15) shows the evolution of a perturbation is along the four basis vectors. For an 
arbitrarily oriented initial perturbation the growth may not be dominated in the direction of 
greatest expansion for short integration times. This can be overcome by sampling multiple 
perturbations in the different directions. A reference point and its neighbors are identified 
and after a finite time their positions in configuration space are computed. The state tran- 
sition matrix can then be computed at each point in the xy plane, by using a central finite 
difference method. For initial perturbations restricted to W x W y subspace, this gives 



,k+l,l (to+T)— Xjjje-iA (tp+T) tti,j,fc,i+i (to+Tj-gj, j,k,i-i (tp+T) 



\rW 



AW x (t ) AW y (t ) 
yi,j,k+i,i{to+T)-yij t k-i,i(to+T) yi,j,k,i+i(to+T)-y i:j:k]l _ 1 (to+T) 
AW* (to) AW y (t ) 



(23) 



The relative velocity sensitivity field , cr(W x , W y ) is given by, 



cr(W x , W y ) = In ^\ ma x{.4>* rW <P,rw) (24) 

Ridges on this sensitivity surface are one dimensional structures similar to LCS. The ridges in 
the maximal sensitivity field o~(W x , W y ), partition the relative velocity subspace. We applied 
the above procedure to a cellular flow. 

A note on the terminology - The field measuring the sensitivity of the final location of 
particles in configuration space with respect to perturbations in initial relative velocity is 
analogous to the FTLE field, but not identical. To obtain the true FTLE field, one would 



have to compute the 4x4 state transition matrix, Using the notation of eq (23) 



$ = Y T Y w (25) 

1 %r <PWW ' 



The FTLE field is then given by a(x, y, W x , W y ) — t4- In yj\ ma x ($*$)• Ridges in this field 
are 3 dimensional structures and represent the true LCS. 

4 Example: cellular flow 

This flow is described by the stream function 

ip(xj y,t) — a cos x cos y (26) 

The velocity field is given by, 

u = — a cos x sin y (27) 
v = a sin x cosy (28) 
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There are heteroclinic connections from the stable and unstable manifolds of the fixed 
points (2n + 1)|, shown by the arrows in Figure [TJ which are also the boundaries of the 
cells. These coincide with LCS/, the LCS for the fluid velocity field. The LCS/ is to be 
distinguished from the LCS of the inertial particle in the full four dimensional phase space. 
By choosing initial perturbations of the form given by eq. (21) at different points along 
a streamline, we follow how these perturbations grow in the xy plane by integrating the 
particle trajectories numerically from which the sensitivity field is computed. Figure [2] shows 
the sensitivity field computed for initial perturbations in the relative velocity subspace, at 
different points on the streamline ip(x, y, t) = 0. The ridges in this field have high values of 
sensitivity. It can be seen that there is a continuous variation in the ridges of the sensitivity 
field with respect to the initial (x, y) coordinates. In each case the sensitivity field at a given 
point depends on the underlying LCS / of the fluid flow. 

The ridges in the sensitivity field have meaningful information about the dynamics of 
inertial particles even when computed at points far from the saddle points of the fluid flow. 
This is shown in Figure [3|a) which is the sensitivity field computed at (x,y) = ^). The 
ridges in the sensitivity field partition the relative velocity subspace according to the final 
location of particles. In Figure |3^b) the ridges in the partial FTLE field are used to identify 
regions in the relative velocity subspace, that produce qualitatively different trajectories. 
Particles that start at the same physical location, but are in different regions of the relative 
velocity subspace, are neatly separated from particles that started in other regions, as shown 
in Figure [31(c) . Thus the ridges in the sensitivity field, have the property of a separatrix. 



5 Segregation of particles by Stokes number 

Eq. ^ can be diagonalised as 

dW d ( -\-n 
dt ~ V A ~ 

where A are the eigenvalues of the Jacobian of the fluid velocity field. If \i = fS't" 1 , is very 
large, then both the components of Ad would decay. For low values of /x, one component of Ad 
would grow. Therefore the dynamics of an inertial particles depend on the value of /i, that is 
on the Stokes number. It is reasonable to expect that the computations of the sensitivity of 
the particles location to the initial relative velocity also would depend on the Stokes number. 
That this is indeed the case is shown by the computations of the sensitivity field for a particle 
with Stokes number 0.1 for the time independent flow, as shown in Figure |4j(a) . The thick 
lines are the ridges in sensitivity field for particles with St = 0.1 and the hatched lines are 
those of St = 0.2. It can be seen that though the structure of the sensitivity field field is 
similar, the ridges are present at different locations in the relative velocity subspace. This 
fact can be exploited to design a process to separate particles by their Stokes number. In 
this section we illustrate a simple procedure for doing this. 

The ridges of the sensitivity fields computed for the two different particles of Stokes 
number 0.1 and 0.2 respectively are overlain in the same plot, as shown in Figure |4j The 
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subdomain of the relative velocity subspace sandwiched between the ridges of the sensitivity 
fields of the two types of particles form a zone of segregation. One such sample zone is shown 
in grey in Figure[4j Two particles with St = 0.1 and St = 0.2 with initial coordinates (x, y) = 



^) and the same initial relative velocity, belonging to this region, have trajectories that 
separate in the physical space. To illustrate this, the trajectories of five hundred particles 
of each Stokes number, starting at the same initial physical point (x, y) — ^) and with 
initial relative velocities values belonging to the grey region were computed. The position 
of these particles is plotted as a function of time to show that the particles are completely 
segregated into two different cells. 

The above procedure can be applied to other regions sandwiched between the ridges of 
the partial FTLE of the two different types of particles. 

6 Robustness of the sensitivity field to perturbations 
in the stream function 



The time- independent flow given the stream function in eq. (26) is perturbed by making it 



weakly time dependent. The modified fluid flow is given by the stream function 

i[)(x,y, t) = a cos (x + bsmuut) cosy (30) 
The velocity field is given by, 

u = — a cos {x -\- b sin out) sin y (31) 



v = asm (x + b sin cut) cosy (32) 

For time dependent systems the location of the LCS depends on the choice of initial 
time. For the computation of the sensitivity field, the location of ridges in the relative 
velocity subspace depend on the initial spatial coordinates of the particle as well as the 
initial time. However our computations show that the dependence of the ridge structure on 
the initial time is weak. Figure [5] shows the ridges in the sensitivity field. As the initial time 
is increased, it is seen that there is a 'squeezing' of the sensitivity field in some regions of the 
relative velocity subspace. A comparison with Figure [4] and Figure [3] shows that the ridge 
locations in the sensitivity field remain qualitatively the same, for the three cases Figure [5] 
where the initial time is small. This offers a numerical evidence that the sensitivity field is 
robust to small perturbations in the fluid velocity. 



7 Conclusion 

The dynamics of inertial particles in a fluid flow can exhibit sensitivity to initial conditions. 
The finite time Lyapunov exponent can be used to characterize this sensitivity. The LCS 
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obtained from the ridges of the FTLE field offers a systematic method to identify qualitatively 
different regions of the phase space. We demonstrated that even a reduced one dimensional 
ridge in the sensitivity field contains important information about the sensitivity of the 
spatial location of particles to initial relative velocity. The Stokes number, and by implication 
the size of the particle, is an important parameter that governs the clustering behavior of 
particles for a given flow. This property can be exploited to make particles of different sizes 
cluster in different regions of the fluid and thus separate them. For the more general case of 
non neutrally buoyant particles, the density of the particles could play a similar governing 
role as the Stokes number. One could therefore design flows that can fractionally separate 
particles for a range of inertial parameters. 
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Figure 1: Streamlines of ip = a cos x cosy form an array of cells. The arrows indicate the heteroclinic 
trajectories connecting the fixed points of the velocity field formed by ip. For this velocity field, the heteroclinic 
trajectories coincide with the LCS 
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Figure 2: Ridges in the sensitivity field for tp = a cos x cosy. Initial spatial position varies from (x ,y ) = 
(n/2, n/2) to (0, 7r/2), at the points shown in (1), along tp = 0, at intervals of 0.057T. The plots show a smooth 
variation in the structure of the ridges in the sensitivity field. Parameters , a = 100, St = 0.2, integration 
time T = 0.24 
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(a) 




Fi guxe 3i (a) Ridges in the sensitivity field, partition the velocity subspace into regions of distinct qualitative 
dynamics. Three such partitioned regions are shown, (b) Particles starting with relative initial velocities 
belonging to distinct partitions in the relative velocity subspace are segregated into different cells in the xy 
plane. Time of integration T = 0.24, St = 0.2. The initial position of all particles is (xojZ/o) = (37r/8, 37r/8), 
shown by the x marker. 
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(b)T = 0.005 (c)T = 0.030 (d)T = 0.060 
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(e)T = 0.085 (f)T = 0.110 (g)T = 0.135 




(h)T = 0.160 (i)T = 0.185 (j)T = 0.210 



Figure 4: (a) Ridges in the sensitivity field for particles with St =0.2 (hatched) and St=0.1(thick) re- 
spectively. The initial position of all particles is (x n , y ) — (37r/8, 37r/8). The grey patch is a sample region 
sandwiched between the ridges of the two Stokes numbers, (b) Particles starting at (x ,y ) = (3tt/8,3tt/8) 
and relative initial velocity belonging to the grey patch, are separated into different cells in the xy plane. 
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Figure 5: (a) Ridges in the sensitivity field for the stream function ip(x,y,t) = a cos (x + b sin cot) cos y for 
(xoiUo) = (37r/8, 37r/8). The hatched lines and the thick lines are the ridges corresponding to St = 0.2 and 
St = 0.1 respectively. Parameters, a = 100, b = 0.25, ui — 1, Integration time = 0.24. Initial time T = 0. (b) 
T = 0.25 (c) T = 0.5 
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